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ABSTRACT 

We have carried out 3-D numerical simulations of the dynamical bar 
instability in a rotating star and the resulting gravitational radiation using 
both an Eulerian code written in cylindrical coordinates and a smooth particle 
hydrodynamics (SPH) code. The star is modeled initially as a polytrope with 
index n = 3/2 and T rot /|VK| « 0.30, where T rot is the rotational kinetic energy 
and \W\ is the gravitational potential energy. In both codes the gravitational 
field is purely Newtonian, and the gravitational radiation is calculated in the 
quadrupole approximation. 

We have run 3 simulations with the Eulerian code, varying the number of 
angular zones and the treatment of the boundary between the star and the 
vacuum. Using the SPH code we did 7 runs, varying the number of particles, the 
artificial viscosity, and the type of initial model. We compare the growth rate 
and rotation speed of the bar, the mass and angular momentum distributions, 
and the gravitational radiation quantities. We highlight the successes and 
difficulties of both methods, and make suggestions for future improvements. 

Subject headings: hydrodynamics — methods: numerical - instabilities - 
radiation mechanisms: gravitational 



1. Introduction 

Many of the most interesting astrophysical systems can be described by the equations 
of hydrodynamics coupled to gravity. As computers grow more powerful, new numerical 
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techniques are being developed and computer simulations of astrophysical systems are 
gaining importance. In fact, numerical modeling may be the only means of getting detailed 
understanding about certain phenomena. Each numerical method has its own strengths and 
weaknesses, and the choice of the most suitable method depends in part on the physical 
system being studied. Therefore, it is important to understand the behavior of the various 
techniques in different situations. 

One area in which numerical simulations play a key role is the modeling of astrophysical 
sources of gravitational radiation. With the prospect of several gravitational wave detectors 
becoming operational within a decade (e.g. Abramovici, et al. 1992; Bradaschia, et al. 1990), 
the detailed modeling of these sources has a high priority. For example, global rotational 
instabilites that arise in collapsing or compact stars can potentially produce detectable 
amounts of gravitational radiation. A rapidly rotating stellar core that has exhausted 
its nuclear fuel and is prevented from collapsing to neutron star size by centrifugal forces 
could become unstable and possibly shed enough angular momentum to allow collapse to a 
supernova (Thorne 1995). Also, a neutron star that is spun up by accretion of mass from a 
binary companion could potentially reach fast enough rotation rates to go unstable (Schutz 
1989). Since these sources are time-dependent, nonlinear, and fully 3-dimensional systems, 
calculating the gravitational radiation they produce requires numerical simulations. 

Global rotational instabilities in fluids arise from nonradial "toroidal" modes e ±trn<p , 
where tp is the azimuthal coordinate and m = 2 is known as the "bar mode" . They can be 
parametrized by 



where T rot is the rotational kinetic energy and W is the gravitational potential energy 
(Tassoul 1978; Shapiro k Teukolsky 1983; Durisen & Tohline 1985; Schutz 1986). We 
concentrate on the bar instability since it is expected to be the fastest growing mode. This 
instability can develop by two different physical mechanisms. The dynamical bar instability 
is driven by Newtonian hydrodynamics and gravity. It occurs for fairly large values of 
the stability parameter r > r d and develops on a timescale of approximately one rotation 
period. The secular instability arises from dissipative processes such as gravitational 
radiation reaction and operates in the range r s < r < t&. It develops on a timescale of 
several rotation periods or longer (Schutz 1989). The constant density, incompressible, 
uniformly rotating Maclaurin spheroids have r s m 0.14 and m 0.27. For differentially 
rotating fluids with a polytropic equation of state, 



where n is the polytropic index and K is a constant that depends on the entropy, early 
studies indicated that the secular and dynamical bar instabilities should occur at about 





P = Kp T = Kp 1+1 / n , 



(2) 
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these same values of r (Shapiro & Teukolsky 1983; Durisen & Tohline 1985; Managan 
1985; Imamura, Friedman, & Durisen 1985). Recent work by Imamura, et al. (1995) shows 
that both the angular momentum distribution and, to a lesser degree, the polytropic index 
affect the value of r at which the m = 2 secular instability sets in. For the dynamical 
bar instability Pickett, Durisen, & Davis (1996; hereafter PDD) demonstrate that, for 
n = 3/2 polytropes, the m = 2 dynamical stability limit r d pa 0.27 is valid for initial angular 
momentum distributions that are centrally condensed and similar to those of Maclaurin 
spheroids. However, for angular momentum distributions that produce somewhat extended 
disk-like regions, both one- and two-armed spiral instablities appear at considerably lower 
values of r. 

As a first step toward understanding realistic sources we are simulating the gravitational 
radiation emitted when a rapidly rotating star, modeled initially as a polytrope with 
n = 3/2 (r = 5/3), becomes dynamically unstable. Newtonian gravity is used, and the 
gravitational radiation produced is calculated in the quadrupole approximation. The back 
reaction of the gravitational radiation on the star is not included. We have chosen the 
case r pa 0.30, which is just above the dynamical stability limit and so might reasonably 
approximate a star that spins up (due to collapse or accretion) and goes unstable. This case 
has also been studied numerically and analyzed using the linearized tensor virial equations 
(TVE; see Chandrasekhar 1969) by Tohline, Durisen, & McCollough (1985; hereafter 
TDM), so their results are available for comparison. 

We have carried out simulations of the dynamical bar instability using two very 
different computer codes, each based on numerical techniques actively used in astrophysics. 
One of these is a 3-D Eulerian hydrodynamics code written in cylindrical coordinates 
with monotonic advection. The other is a smooth particle hydrodynamics (SPH) code 
with variable smoothing lengths and individual particle timesteps. Since the SPH code is 
Lagrangian, gridless, and fully adaptive, it is intrinsically very different from the Eulerian 
code. By running the same calculation on these two codes, we hope to gain a better 
understanding of the relative merits of these methods in modeling the dynamical bar 
instability. 

An earlier comparison between the results of using Eulerian and SPH codes to model 
the dynamical instability was carried out by Durisen et al. (1986; hereafter DGTB). They 
used rapidly rotating polytropes with n = 3/2 and considered the cases r pa 0.33 and 
t pa 0.38. Since they were studying this instability in the context of star formation, they did 
not calculate the gravitational radiation generated. They used two different Eulerian codes, 
one with cylindrical coordinates (the same one used by TDM) and the other with spherical 
coordinates. Both of these used the diffusive donor cell advection and fairly low resolution. 
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The SPH code used a smoothing length that was the same for all particles and varied in 
time, and a fairly small number of particles. Our study takes advantage of more modern 
and accurate numerical methods, and focuses on the gravitational radiation generated by 
the dynamical instability in compact stars. 

This paper is organized as follows. In § ^ we briefly describe the numerical techniques 
used in the two codes, and in § [T] we discuss the calculation of gravitational radiation using 
the quadrupole approximation. The initial conditions are presented in § [TJ. The results of 
modeling the bar instability using the Eulerian code are given in § [5], and the results of 
using the SPH code in § We compare the Eulerian and SPH results in § |7] and present 
our conclusions in § ^ 



2. Numerical Techniques 

Both of the computer codes used in this study solve the equations of hydrodynamics 
coupled to Newtonian gravity. The matter is taken to be a perfect fluid with equation of 
state 

P = (T - l)pe, (3) 

where e is the specific internal energy. Each code has been subjected to a variety of tests to 
insure its accuracy and stability. In this section we present a brief description of each code, 
referring the reader to the literature for further details. 



2.1. Eulerian Code 

We use the 3-D Eulerian hydrodynamics code developed by Smith, Centrella, & 
Clancy (1994; see also Smith 1993; Clancy 1989). This code is written in cylindrical 
coordinates (zcr, z, (p) with variable spatial zoning. This is useful for representing rotating 
configurations, including bars, toroids, and more complicated geometries, all of which may 
exhibit substantial rotational flattening. The hydrodynamical equations are solved using 
time explicit differencing with operator splitting (Wilson 1979; Bowers & Wilson 1991). 
Although the code has the option of allowing the grid to move in the w and z directions, 
for simplicity we hold both grids fixed for the models presented in this paper. We impose 
reflection symmetry through the equatorial plane z = 0, and calculate the full range of the 
angular coordinate <p : — 2n. 
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In Eulerian hydrodynamics fluid is transported from one grid zone to another, and it is 
important to obtain an accurate value for the quantity crossing the zone face. The simplest 
such advection scheme is the donor cell method, which is only accurate to first order and 
produces large numerical diffusion (Bowers & Wilson 1991). To achieve better accuracy 
and less diffusion, the advection terms can be updated using an interpolation method 
that preserves monotonicity in the quantity being advected. This code uses a monotonic 
advection scheme developed by LeBlanc (Clancy 1989; Bowers & Wilson 1991), with all 
spatial finite differences in the advection phase being second order. The spatial differencing 
in the Lagrangian phase is first order except for the "PdV" term, which is second order. 
The code uses first order differencing in time with operator splitting; in general, this results 
in a scheme that is somewhat better than first order in time, but not quite second order. 
The method of consistent advection is used for the angular momentum transport (Norman, 
Wilson, & Barton 1980; Norman & Winkler 1986). Shocks are handled using a standard 
artificial viscosity. For the bar instability runs presented in this paper, shocks occur during 
the later stages of the evolution, when the spiral arms expand supersonically and merge. 
This shock heating and dissipation generates entropy. 

The Newtonian gravitational potential is calculated by solving Poisson's equation on 
the cylindrical grid, with the boundary conditions at the edge of the grid specified using a 
spherical multipole expansion. In finite difference form this becomes a large, sparse, banded 
matrix equation which we solve using a preconditioned conjugate gradient method with 
diagonal scaling (Press, et al. 1992; Meijerink & Van Der Vorst 1981). This is a simple and 
efficient method that requires a minimum of memory overhead, since it does not need to 
store the entire matrix being inverted and takes advantage of existing arrays already set 
aside for temporary storage in the code. Comparison tests with other sparse matrix solvers 
showed that this method produces solutions with the same accuracy using significantly less 
CPU time (Smith, Centrella, & Clancy 1994). Such memory and time considerations are 
very important for the successful implementation of a fully 3-D Eulerian code. 



2.2. TREESPH 



SPH is a gridless Lagrangian hydrodynamics method that models the fluid as a 
collection of fluid elements of finite extent described by a smoothing kernel (Lucy 1977; 
Gingold & Monaghan 1977; see Monaghan (1992) for a review). We have used the 
implementation of SPH by Hernquist & Katz (1989) known as TREESPH. In this code 
each particle is assigned a smoothing length which is allowed to vary in both space and 
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time, thereby achieving roughly the same level of accuracy in all regions of the fluid. The 
use of these variable smoothing lengths as well as individual particle timesteps makes the 
program adaptive in both space and time. TREESPH has the option of evolving either the 
thermal energy or a function of the entropy. Hernquist (1993) has shown that for adaptive 
SPH, in which smoothing lengths vary in time, certain types of errors do not show up in 
the total energy if the thermal energy equation is evolved. However, if the entropy equation 
is used, then conservation of total energy is a good indicator of the global accuracy of the 
calculation. We have chosen to evolve the entropy equation. 

The gravitational forces in TREESPH are calculated using a hierarchical tree method 
(Barnes & Hut 1986) optimized for vector computers (Hernquist 1987). The particles are 
first organized into a nested hierarchy of cells, and the mass multipole moments of each 
cell up to a fixed order, usually quadrupole, are calculated. In computing the gravitational 
acceleration of a particle, it is allowed to interact with different levels of the hierarchy in 
different ways. The force due to neighboring particles is computed by directly summing the 
two-body interactions. The influence of more distant particles is accounted for by including 
the multipole expansions of the cells which satisfy the accuracy criterion at the location 
of each particle. In general, the number of terms in the multipole expansions is small 
compared to the number of particles in the corresponding cells. This leads to a significant 
gain in efficiency, and allows the use of larger numbers of particles than would be possible 
with methods that simply sum over all possible pairs of particles. 

As a Lagrangian method SPH is attractive because the computational resources can be 
concentrated where the mass is located, rather than spread over a grid that can be mostly 
empty. In addition the numerical algorithms are simpler and, in general, considerably 
easier to implement than standard Eulerian methods. SPH has been applied to a variety 
of astrophysical problems and is gaining in popularity. However, it is still a relatively 
new method, and less is known about its behavior in various situations than the Eulerian 
methods, which have been developed and used by a much larger number of researchers over 
the decades. Comparison studies such as this one are therefore of considerable interest. 



3. Calculation of Gravitational Radiation 

We calculate the gravitational radiation produced in these models using the quadrupole 
approximation, which is valid for nearly Newtonian sources (Misner, Thorne, & Wheeler 
1973). Since the gravitational field in both codes is purely Newtonian, we calculate only the 
production of gravitational radiation and do not include the effects of radiation reaction. 
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The spacetime metric can be written 

sv = v^u + V> ( 4 ) 

where /i, v — 0, 1, 2, 3, r\^ v = diag(— 1, 1, 1, 1) is the metric of flat spacetime, and \h^ u \ <C 1. 
The gravitational waveforms are given by the transverse-traceless (TT) components of the 
metric perturbation hij, 

r/ i J T = 2 4 TT , (5) 

where 

% = J P ( x i x j ~ l S ij r ' 2 ) rf3r ( 6 ) 

is the trace-reduced quadrupole moment of the source. Note that we use units in which 
c — G — 1 only when discussing the gravitational radiation quantities. Here r 2 = x 2 + y 2 + z 2 
is the distance to the source, spatial indices i,j = 1,2,3, and a dot indicates a time 
derivative d/dt. For an observer located on the axis at 9 = 0, tp = in spherical coordinates 
centered on the source, the waveforms for the two polarization states take the simple form 

^■+,axis ~\ixx ^yy)i C' 7 ) 

2 •■ 

^x,axis ~Ixy (8) 

The gravitational wave luminosity is defined by 
and the angular momentum lost through gravitational radiation is 



dJi _ 2 
~dT ~ 5 

where the superscript (3) indicates 3 time derivatives, there is an implied sum on repeated 
indices, and the angle-brackets indicate an average over several wave periods. For these 
burst sources such averaging is not well-defined; therefore we display the unaveraged 
quantities ^if^lf^ and \^ijk^jm^km De l° w - The energy emitted as gravitational radiation is 

AE = J Ldt (11) 

and the angular momentum carried away by the waves is 

AX, = J (d.Ji/dt)dt. (12) 
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The expressions given above for calculating gravitational radiation in the quadrupole 
approximation are all functions of at least the second time derivative of f^. The standard 
quadrupole formula consists of using the definition (|5|) for in these equations. In an 
Eulerian code, iy may be calculated directly by summing over the grid and the time 
derivatives may be taken numerically. However, this successive application of numerical 
time derivatives can introduce a great deal of noise into the calculated quantities, especially 
when the time step varies from cycle to cycle. 

To reduce this problem, Finn & Evans (1990) have developed two partially integrated 
versions of the standard quadrupole formula that eliminate one of the time derivatives; 
they call these the momentum divergence and the first moment of momentum formulae. 
In these expressions, fjj is calculated directly by integrating fluid quantities over the grid. 
Since eliminating one numerical time derivative in a finite difference code greatly increases 
the signal-to-noise ratio, both of these methods significantly reduce the high frequency 
numerical noise and produce much cleaner waveforms than the standard quadrupole 
formula. 

Both of these formulae from Finn & Evans (1990) have been implemented in the 
Eulerian hydrodynamics code to calculate the gravitational radiation quantities; details are 
given in Smith, Centrella, & Clancy (1994). Since the resulting waveforms are very similar, 
we only show the waveforms obtained using the first moment of momentum expression in 
this paper. Note that this expression gives f^; the waveform requires taking another time 
derivative to obtain f^, and the luminosity requires still another derivative. When these 
derivatives are calculated numerically the signals can still be dominated by noise, especially 
when the time step is changing significantly from cycle to cycle. This problem was solved 
by passing the data through a filter to smooth it after was calculated, and again after 
each numerical derivative was taken. These techniques produce smooth profiles for both 
the waveforms and luminosities, as shown below. 

The gravitational radiation is computed in TREESPH using the method of Centrella & 
McMillan (1993) to calculate analytically from the SPH equations of motion. With this, 
the gravitational waveforms are calculated directly using the particle positions, velocities, 
and accelerations which are already available in the code. The resulting waveforms are very 
smooth functions of time and require no filtering or smoothing to remove numerical noise. 
However, the luminosity and angular momentum lost by gravitational radiation do contain 
the time derivative of the particle acceleration, which is taken numerically and therefore 
introduces some noise. We have chosen to smooth the luminosity data for the SPH runs 
using simple averaging over a fixed interval of 0.1£d centered on each point. Here, is the 
dynamical time defined in equation (|14]) below. In general this procedure makes a negligible 
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change in the integrated luminosity (p~T|) , which gives the energy emitted as gravitational 
radiation, and produces very smooth profiles (Centrella & McMillan 1993). The profiles of 
the angular momentum lost to gravitational radiation are not smoothed. 



4. Initial Models 

The initial conditions for our simulations are rotating axisymmetric equilibrium models 
having r ~ 0.30 and polytropic index n = 3/2. The bar instability then grows from 
nonaxisymmetric perturbations of these equilibrium spheroids. In this section, we briefly 
describe the construction of these initial models and their representations in the Eulerian 
code and in TREESPH. 

The self-consistent field method of Smith Sz Centrella (1992) is used to generate the 
axisymmetric equilibrium models. This is based on the earlier work of Ostriker & Mark 
(1968) and Hachisu (1986), and derives from an integral formulation of the equations of 
hydrodynamic equilibrium which automatically incorporates the boundary conditions. We 
use a cylindrical grid (w, z) with uniform zoning. An initial "guess" density distribution is 
given, and the gravitational potential is calculated using a Legendre polynomial expansion to 
solve Poisson's equation (Hachisu 1986). A rotation law of the general form j(m) = j(m(zu)) 
is specified, where j(m) is the specific angular momentum and m(w) is the mass interior to 
the cylinder of radius vo (Ostriker & Mark 1968). Following the convention of earlier work 
(Bodenheimer & Ostriker 1973; TDM; DGTB; Williams & Tohline 1987, 1988) we use the 
rotation law for the uniformly rotating, constant density Maclaurin spheroids, which can be 
written in the dimensionless form 

h(m) = —j(m) = 1(1- (1-m) 2 / 3 ), (13) 

J 

where J is the total angular momentum and M is the total mass. Since polytropes do 
not have constant density, this produces differentially rotating models. The rotation law is 
used to calculate a rotational potential, which is then used with the gravitational potential 
to compute an improved density distribution. This process is iterated until convergence is 
achieved. 

The freely specifiable quantities in this method are the dimensionless rotation law 
h(m), the polytropic index n, and the axis ratio R p /R eq , where R p is the polar radius and 
R eq is the equatorial radius of the initial model. To get a dimensional model, we also specify 
the maximum density and the entropy, which is given by the constant K in the polytropic 
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equation of state (0). Upon convergence to a solution of the equations of hydrodynamic 
equilibrium, this procedure gives the density p(w), the angular velocity Q(zu), the total 
mass M, the total angular momentum J, and the stability parameter r. Since r is not 
specified initially, some experimentation with the axis ratio is generally necessary to achieve 
a desired value of r. In constructing our model with r ~ 0.30 we were guided in our choice 
of input parameters by the values given in TDM. We found that using R p /R eOL = 0.205 
gives t = 0.301; note that this configuration is highly flattened due to rotation. The central 
rotation period for this model is 2.15£d and the rotation period for a point on the equator 
is 6.90£d, where 

{ R 3 \ 1/2 

is the dynamical time for a sphere of radius R eq . To construct the initial model, we used 
a uniform cylindrical grid of N m = 65 radial zones and N z = 23 axial zones. The mass 
distribution extended out to zone 61 in the w direction and to zone 19 in the z direction. 
This model required 50 iterations to converge to a solution with a tolerance of 10~ 10 and 
used 400 seconds of CPU time on the Cray C90 at the Pittsburgh Supercomputing Center 
(PSC). The accuracy of this initial equilibrium model can be measured using the virial 
relation. Let 

VC= ^±M , (15) 

where II = / PdV is the volume integral of the pressure (Hachisu 1986). For our initial 
model, VC = 7.4 x 10" 4 . 

To evolve this model with the Eulerian code, we first interpolate it onto the non-uniform 
grid used in that code. We trigger the bar instability by imposing a random perturbation 
with amplitude 10 -3 on the density in each grid zone (cf. TDM). The instability then grows 
from this relatively low noise level, with the start of the gravitational wave burst occurring 
at ~ 15£d- 

We have used two different methods to convert the density p(tu, z) and angular velocity 
VL{w) produced by the self-consistent field method into a particle model to be evolved with 
TREESPH. Both methods use equal-mass particles. 

The first technique is a simple random or "rejection" method (Press, et al. 1992; 
Centrella & McMillan 1993) that randomly distributes particles within the probability 
distribution p(w,z), and then assigns the appropriate angular velocity. Since the particles 
are accepted into the model randomly, and thus independently of each other, this method 
results in both positive and negative density fluctuations about the target p(w,z). These 
fluctuations are relatively large. They trigger the bar instability, with the gravitational 
wave burst starting immediately (see model T6 below). 
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Since the perturbations imposed on the Eulerian models at the initial time are 
significantly smaller than those produced by the random particle method, we developed a 
technique for producing a quieter, "cold" initial particle model. In this method a set of 
equipotential surfaces is calculated for the equilibrium model produced by the self-consistent 
field method using the gravitational and rotational potentials. The mass distribution for 
this model is specified by calculating the mass interior to these equipotential surfaces. To 
create a particle representation, we start by placing particles within the surface boundary 
of the star at uniform Cartesian coordinates. The mass interior to the equipotential 
surfaces is then used to determine how to relocate these particles to produce the desired 
mass distribution. Particle velocities are assigned using Q(zu). The resulting models have 
considerably less noise, with the gravitational wave bursts beginning at ~ 10£d (Houser & 
Centrella 1995). 



We ran three Eulerian models, labeled El - E3; in all these models, we use N m = 64 
zones in the w direction and N z = 32 zones in the z direction. To maximize both resolution 
and efficiency we use a finer grid in the region initially occupied by the matter and a 
coarser grid outside. The w grid is uniform up to the zone j = 30, and the z grid is 
uniform up to k = 16. The zoning is chosen such that the center of the radial zone 
j = 25 is at the equatorial radius of the initial model R eq and the center of the axial zone 
k = 9 is at the polar radius R p . Outside of this uniformly zoned region, the zone size 
increases linearly with zoning ratios AtOj + i/AtUj = 1.03 and Az^+i/ Az^ = 1.1 The angular 
grid is uniform and covers the range <p : — 2n. For simplicity, the grids are held fixed 
throughout the runs. The grid boundaries for the Eulerian runs are set at w = 3.85R eq and 
z = 1.60-Rcq = 7.91-Rp. This large amount of initially empty space is necessary to provide 
room for the star to expand as the bar mode grows, and to specify the boundary conditions 
for the solution of Poisson's equation accurately (Smith, Centrella, & Clancy 1994). We 
varied the number of angular zones and the vacuum boundary conditions to test the 
effects of these parameters on the bar mode instability. The properties of the Eulerian 
models are summarized in Table [I]. 

To study the development of the bar mode quantitatively, we analyze the density in a 
ring of fixed w and z using a complex Fourier integral 

1 r 27r 



5. Evolution of the Bar Instability: Eulerian Runs 




(16) 
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where m = 2 (TDM). The normalized bar mode amplitude is 

\C\ = \C 2 \/C , (17) 

where C (w, z) = p(zu, z) is the mean density in the ring. The phase angle m is defined by 

(f) m (tu, z) = tan -1 Im(C m )/Re(C m ). (18) 

The phase information can be used to describe global nonaxisymmetric structure 
propagating in the ^-direction. When such a global mode develops out of the initial noise 
we can write 

4>m = V m t, (19) 

where the pattern speed of the mth structure is (cf. Williams & Tohline 1987; PDD) 

W m {u7, z) = — — = — . (20) 
m at m 

Thus, o"2 is twice the bar rotation speed, and the rotation period of the bar is Tb ar = 47r/o" 2 . 

The bar mode amplitude |C| and phase angle 02 have been calculated by TDM using the 
linearized TVE method. This technique gives exact results for small oscillations of uniform 
density, incompressible ellipsoids such as the Maclaurin spheriods (Chandrasekhar 1969). 
It was adapted to study rotating compressible fluids by Tassoul & Ostriker (1968), and 
applied to rotating polytropes by Ostriker & Bodenheimer (1973). For compressible fluids, 
the TVE method gives only approximate results; see TDM for a discussion. Nevertheless, 
it provides a useful point of comparison for the numerical simulations. TDM used the 
Ostriker-Bodenheimer TVE code to calculate the bar mode growth rate and eigenfrequency 
for the case r = 0.301. According to their analysis, the amplitude |C| should grow 
exponentially with time as the instability develops, with <iln \C\/dt = 0.728 ± O.OSSto 1 (we 
have converted from their units). For the eigenfrequency they find a 2 = 1.892 ± 0.094t5 1 - 
The errors quoted by TDM are of the order ±5% and only account for the expected 
inaccuracies in the equilibrium models. 



5.1. Results from our Standard Eulerian Model 



Our standard Eulerian model is E2, which uses N v = 64. Density contours showing 
the development of the bar instability in this model are presented in Figure [1]. The growth 
of the m = 2 mode produces a bar-shaped structure. This rotating bar develops a spiral 
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arm pattern as mass is shed from the ends of the bar. The bar and spiral arms exert 
gravitational torques, causing angular momentum to be transported. The spiral arms 
expand supersonically and merge together, causing shock heating and dissipation in the 
disk surrounding the central core. The system remains highly flattened throughout, and 
evolves toward a nearly axisymmetric final state. 

The distribution of mass m(w) is shown in Figure ^| for the initial time (dot-dashed 
line), the intermediate time t = 21. li^ (dashed line), and the final time t = 34t^ (solid 
line). The angular momentum J(w) is shown for these same three times in Figure |3|. Note 
that J(w) is normalized by the total angular momentum in the system at the time, which 
is less than the initial value due to non-conservation; see Table [I]. We define the core to be 
all matter contained within cylindrical radius w = R eq . In the final state, the core has 96% 
of the mass, and 86% of the angular momentum; see Table |2|. 

The growth of the bar mode for model E2 is shown quantitatively in Figure |], where 
In \C\ is plotted versus time for the ring w = 0.362i? eq in zone j = 10 in the equatorial plane 
z — 0. We also checked the growth of In \C\ at several other values of w and found that the 
growth rate din \ C\/dt is essentially independent of cylindrical radius within the core. This 
shows that the bar mode grows at a well-defined rate (TDM). To determine the bar growth 
rate we fit a straight line through the data points in Figure |] in the time interval during 
which the function In \ C\ is growing linearly with time. The endpoints of this time interval 
are chosen "by eye" . Then, using the definition that a line segment consists of at least 10 
successive points, the slope is calculated by linear regression for all possible line segments 
in this time interval. The average of these slopes is used to determine the growth rate. For 
model E2 we find d\n\C\/dt = 0.58^ . 

To determine the eigenfrequency 02, we plot cos 02 as a function of time, and use a 
trigonometric fitting routine to calculate 02- The function cos 02 is used for simplicity, since 
02 itself is multi- valued due to the tan -1 in equation ([18]). The fit is performed over the 
same interval used to calculate the bar mode growth rate. As a check on this procedure, 
we also use an FFT to calculate the frequency spectrum. For model E2 we obtain the 
eigenfrequency <j 2 = 1.8^ , which gives a bar rotation period T bar ~ 7t D . Comparison with 
Figure |] confirms that the initial exponential growth of the bar mode takes place over 
approximately one bar rotation period. 

We calculate the growth of other Fourier components of the density in this same ring 
using equation ([16]) and normalizing the resulting amplitudes by Co- Figure [5] shows the 
growth of the amplitudes of the components (a) m = 1, (b) m = 3, and (c) m = 4. In each 
plot, the amplitude of the bar mode is shown as a solid line for comparison. As expected, 
the growth of the bar mode dominates the initial stage of the evolution, with the other 
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components becoming important at later times. In particular, Figure |5] (c) shows that the 
m = 4 mode starts growing after the bar mode is well into its exponential growth regime. 
The m = 4 mode also grows exponentially, but at a faster rate d\a(\Ci\/\Co\)/dt = 1.1*5 1 
than the bar mode. Both modes reach their peak amplitudes at about the same time, then 
drop to local minima and grow again. The eigenfrequency of the m = 4 mode is o± = 3.4t5\ 
giving a pattern speed for this mode W4 = 0.85*5 ■ Since the pattern speed of the bar mode 
is W 2 = 0.9*5 1 ~ W4, this suggests that the m = 4 mode is a harmonic of the bar mode, 
and not an independent mode. This agrees with the expectation that if the m = 4 mode 
were an independent mode, then it would grow at a slower rate than the bar mode. See 
Williams & Tohline (1987). 

In addition, Figure [5] (a) shows that the m = 1 disturbance grows to nonlinear 
amplitude after the bar mode amplitude has reached its maximum value. Recent work by 
Bonnell (1994; see also Bonnell & Bate 1994) in the context of star formation shows similar 
behavior. PDD also find that both m — 1 and m = 2 modes arise for certain initial angular 
momentum distributions. These issues should be investigated more closely in future work. 

Figure |] shows that the amplitude of the bar mode peaks at * ~ 22*d and then drops 
to a local minimum at * ~ 26*d- It then rises to a local maximum near * ~ 30*d and 
subsequently drops again. These features can also be seen in the density contours given in 
Figure || as follows. Focus on the second highest contour starting in frame (b). This reaches 
a maximum bar-like shape between frames (d) and (e), and then grows more axisymmetric 
until around the time of frame (g), which is near the time at which the bar mode reaches 
its local minimum. The contour again develops a bar-like shape before becoming more 
axisymmetric. The behavior of the stability parameter r = T rot /|W| is shown as a function 
of time by the solid line in Figure || for comparison, the dashed line gives Itotai/I^l- Note 
that r reaches its minimum value when the amplitude of the bar mode peaks at t ~ 22*d- 
It then rises to a local maximum r ~ 0.27 near the time * ~ 26*d when the bar mode 
amplitude is at its local minimum, and falls as the bar mode amplitude grows again. This 
anti-coincidence of the bar amplitude and r results from the fact that the higher amplitude 
bar has a greater moment of inertia, which reduces the rotational kinetic energy. 

The behavior of the gravitational wave quantities is also strongly linked to the 
amplitude of the bar mode. Figure [7] shows the gravitational waveforms (a) rh + and (b) 
rh x for an observer on the axis at 9 — 0, (p — 0. The gravitational waveforms show 
a strong initial burst that peaks around the same time that the bar mode reaches its 
maximum amplitude, * ~ 22*d- This is followed by a weaker secondary burst that peaks 
around * ~ 30*d, corresponding to the secondary maximum of the bar mode amplitude. 
Figure |8] shows (a) the gravitational wave luminosity L, (b) the energy AE/M emitted as 
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gravitational waves, (c) the rate dJ z /dt at which angular momentum is carried away by 
the waves, and (d) the angular momentum A J z / J Q lost to gravitational radiation. The 
luminosity L and dJ z /dt both show a primary peak at t ~ 22t^ and a secondary peak at 
t ~ 30tD, separated by a local minimum at t ~ 26£d- Some of the interesting gravitational 
wave properties are summarized in Table |3|. 



5.2. Results from Eulerian Models with Different Parameter Values 

Model El is the same as E2 except that the resolution in the <p direction is reduced 
by a factor of two, giving N v = 32 angular zones. We chose to change N v because 
we are primarily interested in the growth of nonaxisymmetric modes, and thus the 
angular resolution is expected to be an important parameter. The bar growth rate is 
d\n\C\/dt = 0.53t D \ which is ~ 9% smaller than in E2. Since both these models have 
reasonably long, well-defined linear growth regions for d In \C\/dt we believe that these 
differences are real and not just the result of the method used to compute them. The 
eigenfrequency obtained for model El is = 1.7^ 

TDM (see also Norman, Wilson, & Barton 1980; Williams & Tohline 1987) showed 
that, for the first-order donor cell advection scheme, the difference between the true growth 
rate and the actual growth rate in the code is given by a numerical diffusion term. The size 
of this diffusive term is proportional to the size of the angular zones A(p. They also showed 
that, at least to first order, the eigenfrequency 02 is not affected by this numerical diffusion. 
Our code uses a monotonic advection scheme which is significantly less diffusive than the 
donor cell method (Bowers & Wilson 1991; Hawley, Smarr, & Wilson 1984), and in fact 
the bar mode growth rates we obtain are larger, and hence closer to the TVE values, than 
those found by Tohline and collaborators. Nevertheless, we expect that increasing the size 
of A(p by decreasing N v will also lower the growth rate in our code, and this is the behavior 
that we find in comparing El and E2. The differences in o"2 between El and E2 are about 
a factor of 2 smaller than the differences in the growth rate. 

The bar mode amplitude in run El peaks at t ~ 22£d, then drops off and oscillates 
around a lower value for ~ lO^D; and begins to grow again. Both the mass m(w) and 
angular momentum J(w) distributions in run El are more spread out than in run E2, 
resulting in a core {w < R eq ) with a smaller mass and angular momentum. See Tables [3] 
and 0. 

The gravitational waveforms for El show a strong initial burst with maximum 
amplitude ~ 7% smaller than in E2. This is followed by some additional waves, but they 
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are not cleanly organized into a secondary burst as in E2. The luminosity L shows both a 
primary and a secondary peak. And, although the maximum luminosity in El is only ~ 5% 
smaller than in E2, the ratio of the amplitudes of the primary and secondary peaks in L is 
~ 20 for El and ~ 9 for E2. 

The mass density within an Eulerian grid zone can never be zero, since this leads to 
divisions by zero in the code. Therefore, "vacuum" regions of the grid actually have a very 
small mass density. However, if allowed to evolve unrestricted, these low density zones can 
attain very high velocities and begin to dominate the timestep calculations. To prevent 
this, special provisions must be made to handle these "vacuum" regions (R. Bowers, private 
communication, 1991). We have chosen to place the following restrictions on low density 
zones. For a grid zone in which the density is below a certain threshold value, the velocity 
is set to zero. The density threshold we use to limit the velocity for our standard run E2 
(as well as for El) is 10~ 7 p maXj i, where p max ,i is the maximum density at the initial time. 
Also, if the density in a zone is < lCr 10 p max ,i, the internal energy is set to a fixed value 
that produces consistency between equations (0) and (|3]). Finally, the density itself is set 
to lCr 15 p maX)i if it is less than this value. These conditions lead to some loss of energy 
and momentum as matter flows into these cells. DGTB report a similar loss of angular 
momentum that they attribute to the zeroing of velocities in low density zones. 

To see how the vacuum restrictions affect the evolution of model, we ran a simulation 
with relaxed vacuum restrictions. Model E3 is the same as E2 except that the thresholds 
specifying the vacuum conditions are less restrictive, and the velocity is never set to zero. 
For run E3, the velocity in a zone is set to the sound speed if it exceeds the sound speed 
and the density is below 10~ 14 p maX) ;. The threshold below which the internal energy is set 
to a fixed low value is 10~ 13 p maXji . The density itself is set to 10~ 15 p maXji if it is less than 
this value, as in E2. 

Table [l] shows that the simulation time covered by E3 is 23.9£d, which is considerably 
less than the 34.0£d covered by E2. Run E3 was not continued beyond this point because 
this would have been too expensive in terms of CPU time. Athough less simulation time is 
covered, E3 takes almost as much computer time as E2, with the last 1.5£d of simulation 
time for E3 using 5 CPU hours even with a relaxed Courant condition to allow 20% larger 
timesteps. Because of the high velocities of the matter in the low density zones, the timestep 
continually decreases throughout the run, and we believe that it would take ~ 100 hours of 
CPU time to run the remaining ~ 10£d to complete E3. This demonstrates the need for 
restrictions on the vacuum zones. 

The less restrictive vacuum conditions contribute to better energy and angular 
momentum conservation. By the end of run E3 the model has lost 1.5% of its initial energy, 
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and 3.4% of its initial angular momentum. At this same time in E2, the model has lost 2% 
of its energy and 4.9% of its angular momentum. 

Although E3 was not run long enough to evolve the model completely, sufficient time 
has elapsed for some useful comparisons. The bar mode amplitude has peaked and spiral 
structure has developed, but the model has not yet settled back to the nearly axisymmetric 
final configuration of E2. The gravitational wave amplitudes have also peaked by this time, 
although the initial wave burst is not yet complete. We can thus compare the bar mode 
properties and the peak gravitational wave signals. To provide a reasonable comparison 
of other properties, we evaluate them for E2 at 23.9£d (these entries in Tables [1| — ^] are 
labeled E') and compare them to the final results for E3. 

Examining the bar mode properties of E3 shows that they are very similar to those 
for E2. The growth rate for E3 is only slightly smaller, and the eigenfrequencies are the 
same to within the limit of our measurement accuracy. This is not surprising, since the 
growth of the bar mode occurs before there is significant expansion of the model. Also, 
these properties are measured within the central bulk of the configuration, so they should 
not be strongly affected by the treatment of the vacuum regions. An examination of the 
density contours, however, shows that the system expands significantly more when the 
vacuum conditions are relaxed. Figure |9] shows density contours for (a) E2 and (b) E3 at 
time t = 23.9£d- The inner two contours are nearly identical, but the spiral arms ejected by 
the spinning star extend out to a larger radius than in E2. 

The gravitational wave pulses are qualitatively similar but model E3 shows somewhat 
reduced amplitudes, with the maximum amplitude ~ 9% smaller than that of E2. The 
difference is more pronounced when we examine the luminosity, with the peak luminosity 
of model E3 ~ 19% smaller than that of E2. The peak gravitational radiation amplitudes 
are thus somewhat sensitive to the treatment of the boundary between the fluid and the 
vacuum, which affects the outer, lower density regions of the star. However, this is not 
expected to be a major factor in situations that are not dominated by expansion, such as 
rotating stellar core collapse. 

After this work was completed, we learned that R. Durisen and collaborators have 
carried out similar calculations using an Eulerian code (PDD). They were able to avoid 
problems arising from the time step becoming too small without inhibiting expansion 
by setting the background density to a value between 1CT 10 and 1CT 7 p max ,i and limiting 
all velocities in the background to less than twice the maximum initial sound speed (R. 
Durisen, private communication). We plan to incorporate their suggestions into our future 
simulations. 
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An alternative means of achieving a better treatment of the interface between the 
star and the vacuum might be to use the piecewise parabolic method (PPM; Colella 
& Woodward 1984; Davies, et al. 1993), which is known to be very good at handling 
discontinuities in the flow. It also has a higher resolution for a fixed number of zones than 
the finite difference scheme used here. Of course, PPM is also more expensive in terms of 
CPU usage. A comparison of this model run on a PPM code would be very interesting. 

We can also compare our results with other numerical calculations. Tohline and 
collaborators used a 3-D Eulerian code written in cylindrical coordinates. Their code does 
not solve an energy equation, and therefore has no way to handle self-consistently the 
shocks that form as the spiral arms expand. Instead, they required that the fluid maintain 
the same polytropic equation of state (0), and hence the same entropy, throughout the 
evolution. They used donor cell advection, which is known to be very diffusive. And, their 
code assumes "7r-symmetry" , which means that the flow is calculated only in the angular 
range < (p < n so that only the even mode distortions are modeled. TDM used iV ro = 31, 
N z = 15 and N v = 32, with the model extending out to zone 24 in the ro-direction and 
zone 9 in the z-direction; this is essentially the same as our resolution for E2 and E3. They 
found a bar mode growth rate of din \C\/dt = 0.22^ and an eigenfrequency cr 2 = 2.1^ ■ 
They demonstrated that the substantial deviation from the TVE result for the growth rate 
is due to the large numerical diffusion in their code; see TDM for details. Williams and 
Tohline (1987) modeled the initial development of the bar instability in a polytrope with 
r = 0.31 using the same code with N m = 32, N z = 32, and N v = 64, again employing 
7r-symmetry. This doubling of the number of angular zones increased the bar mode growth 
rate to dln\C\/dt = 0.49t^ ; the eigenfrequency was o~2 = 1-9^ ■ They found that the 
m = 4 mode starts growing exponentially after the bar mode does, and that it grows at 
a faster rate. Also, the m = 4 pattern moves together with the m = 2 pattern, with the 
maxima locked in phase, implying that the m = 4 pattern is a harmonic of m — 2, and not 
a distinct mode. Finally, PDD used a modified version of Tohline's code which is second 
order in all spatial differences, including advection terms, and second order in time. They 
calculated the development of the bar instability in a polytrope with r = 0.304 using N m = 
64, N z = 16, and N v = 64 without 7r-symmetry. They obtained <im \C\/dt = O.SS^ 1 and 
a 2 = 2.lt^ (PDD). 



6. Evolution of the Bar Instability: SPH Runs 



We ran a series of 7 models using TREESPH, labeled Tl - T7. All of these models 
use equal-mass particles. The initial state for each model was produced using the "cold" 
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method described in § |1] except for T6, in which the random method was used. We vary the 
number of particles N and the linear and quadratic artificial viscosity coefficients a and /3, 
respectively. In all cases, the number of neighbors that contribute to the smoothing kernel 
is chosen to be A/"s = 64. The properties of the SPH models are summarized in Table [|. 



6.1. Results from our Standard SPH Model 



Our standard SPH model is T7, which has N = 32, 914 particles, a = 0.25, and 
(3 = 1.0. Figure ITU] shows all the particles in this model projected onto the equatorial plane. 



The system rotates in the counterclockwise direction. We see the development of the bar 
and then the spiral arm pattern as mass is shed from the ends of the rotating bar. The 
gravitational torques exerted by the bar and spiral arms cause angular momentum to be 
transported. The spiral arms expand supersonically and merge, causing shock heating. The 
system then evolves toward a nearly axisymmetric final state. Throughout its evolution the 
system remains flattened. Figure [11] shows the particle positions projected onto the x — z 
plane at the final time t = 35£d- Note that the disk around the central core is somewhat 
puffed up due to shock heating during the expansion of the spiral arms. 



Density contours in the equatorial plane for model T7 are shown in Figure 12, with the 



frames corresponding to the same ones displayed in Figure [10]. To produce these contours 
we first use kernel estimation to interpolate the density of T7 onto the cylindrical grid 
used for the Eulerian model E2. We then calculate contours for the matter located in the 
equatorial plane. The contour levels are the same as those used in Figure [I] for model E2. 

The mass distribution m(vj) is shown in Figure [13] for the initial time (dot-dashed 
line), the intermediate time t = 18.2£d (dashed line), and the final time t = 35£d (solid line) 



for model T7. Figure |14] shows the angular momentum distribution J(w) for these same 
three times. As before, we define the core to be all matter contained within cylindrical 
radius w = R eq . In the final state the core has 90% of the mass and 72% of the angular 
momentum; see Table [5]. 

To study the development of the Fourier components of the density, we first interpolate 
the particle model onto the cylindrical grid every O.liD- We then use the same procedure 
developed for the Eulerian models and analyze the density in the ring w = 0.362i? eq in 
zone j = 10 in the equatorial plane using equations flltf ) - (]20"D. The growth of the bar 



mode is shown quantitatively in Figure [15], where In \ C\ is plotted versus time. Comparison 



of Figure 15 with Figure f| shows that the region in which In \C\ grows linearly with time 
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is shorter and less clearly defined in run T7 than in run E2. Since the endpoints of this 
linear region are chosen "by eye," there is a somewhat greater element of subjectivity in the 
resulting bar mode growth rate for T7; this applies to all the TREESPH runs. For T7 we 
find that the growth rate is din \C\/dt = O.Slt^ 1 and the eigenfrequency is 02 = l-9i^ ■ 



Figure [16] shows the growth of the amplitudes of the components (a) m — 1, (b) m — 3, 
and (c) m = 4. In each plot, the amplitude of the bar mode is shown as a solid line for 
comparison. As we saw in § [5]. the growth of the bar mode dominates the initial stage of 
the evolution. The m = 1 and m = 3 components do not exhibit significant growth. The 
m = 4 mode starts growing after the bar mode is well into its exponential growth regime, 
and grows at a faster exponential rate cZ ln( IC4I / |C | ) — l-ltr> ■ Both the m = 2 and 
m = 4 modes reach their peak amplitudes at about the same time, then drop to local 
minima and grow again. The eigenfrequency of the m = 4 mode is 04 = S.St^ 1 , which gives 
a pattern speed for this mode W4 = 0.95^. Since the pattern speed of the bar mode is 
W2 = 0.95^ = W4, this implies that the m = 4 mode is a harmonic of the bar mode, and 
not an independent mode. 

Figure [15] shows that the amplitude of the bar mode peaks at t ~ 18.5£d and then 
drops to a local minimum at t ~ 25£d- The amplitude then begins to rise sharply again 
and levels off around t ~ 30£d- Compare this behavior with that seen in the contour plots 
in Figure [12] by focussing on the innermost contour starting in frame (b). This shows the 
pronounced initial development of the bar, with the maximum elongation of the contour 



occurring around the time of frame (d), in agreement with Figure [15| This contour then 
becomes more axisymmetric until around the time of frame (g), after which it elongates 
again and then grows more axisymmetric. Thus, in run T7, this contour becomes nearly 
axisymmetric before the bar mode amplitude reaches its local minimum. For comparison, 
the behavior of the stability parameter r is shown as a function of time by the solid line in 
Figure [L?|; the dashed line gives T to tai/|W^|- Note that r reaches a local minimum when the 
amplitude of the bar mode peaks at t ~ 18£d- It then rises to a local maximum r ~ 0.27 at 
t ~ 22£d and drops off again. 

The gravitational radiation quantities also exhibit features corresponding to the 



behavior of these modes. Figure [18] shows the gravitational waveforms (a) rh + and (b) 
rh x for an observer on the axis at 9 = 0, (p = 0; cf. Table [(]. The gravitational waveforms 
show an initial burst that peaks at t ~ lSt-a, corresponding to the initial growth of the 
bar instability. After the peak, the amplitude drops off until t ~ 22t^ and then stays at a 
nearly constant value for about one bar rotation period, during which time the bar mode 
amplitude reaches its local minimum. Recall that Tb ar = 2Tqw, where Tqw is the period 
of the gravitational waves. The wave amplitude then drops again at t ~ 30to and stays 
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at a nearly constant amplitude for about another bar rotation period before the run ends. 
These low amplitude waves are produced by the rotating, slightly non-axisymmetric core 
in the final state; cf. Figure |I^. Figure shows (a) the gravitational wave luminosity L, 
(b) the energy AE/M emitted as gravitational waves, (c) the rate dJ z /dt at which angular 
momentum is carried away by the waves, and (d) the angular momentum AJ z /Jq lost to 
gravitational radiation. The luminosity L shows a broad primary peak centered around 
t ~ 18t-£,, followed by a drop to a local minumum around t ~ 23£d and then a secondary 
feature at t ~ 25to- The signal then drops again to a nearly constant level for t ^ 30tD- 



6.2. Results from SPH Models with Different Parameter Values 

To understand how changing the resolution affects the behavior of the model, compare 
runs Tl, T2, and T3 with the standard run T7. These runs differ only in the total number 
of particles N; see Table In all these cases, the bar and spiral arms develop as in T7. 
Plots of the particles projected onto the equatorial plane appear visually similar, except that 
the extent of the spiral arm pattern increases somewhat as N increases due to the larger 
number of particles available to resolve the outer regions. Table [5] shows that the properties 
of the cores of these runs are all similar. The behavior of the Fourier components of the 
density in these runs is also similar, except that the results for Tl are much noisier due to 
its very low resolution. Although Tl has the largest bar mode growth rate, we attribute 
this to the lack of a clearly defined linear growth region for din \ C\/dt and therefore do not 
consider it to be a reliable indicator of the accuracy of this model. 

The gravitational wave quantities show definite trends with particle number N. As N 
decreases, the amplitude of the burst goes down and the structure of the waveforms after 
the burst becomes less distinct. The peak amplitude of the luminosity L also decreases, 
and the secondary peak becomes a plateau and then disappears into noise for the poorly 
resolved run Tl. For the sequence of models T7, T3, and T2, each run has roughly half the 
number of particles as the previous one. Quantitatively, the amplitude of the burst goes 
down by ~ 10% between model T7 and T3, and another ~ 10% between T3 and T2. The 
peak values of L and dJ z /dt decrease by ~ 30% between T7 and T3, and by another ~ 20% 
between T3 and T2. See Table ^. 

Thus, the larger the number of particles N, the larger the amplitudes of the 
gravitational wave signals. It is clear from Table ^ that the values of these amplitudes 
have not yet converged. One possible means of achieving convergence is simply to increase 
N. However, doubling the number of particles increases the effective grid resolution of the 
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model by only ~ 2 1 / 3 , since we are working in 3-D. Perhaps a better method would be to 
use non-equal-mass particles, with the lower mass particles distributed in the lower density 
regions and thus increasing the resolution in the outer parts of the star (e.g. Monaghan & 
Lattanzio 1985). 

In numerical simulations, artificial viscosity terms are typically added to the 
momentum and thermal energy equations to provide a dissipative mechanism that converts 
the energy jump across the shock into heat (Bowers & Wilson 1991). This smooths out the 
discontinuities that occur in shock fronts while satisfying the Rankine-Hugoniot relations, 
which specify the conservation of mass, momentum, and energy across the shock (Landau 
& Lifshitz 1959). If there is no artificial viscosity, the kinetic energy of the matter passing 
through the shock is not correctly converted into heat and there can be large post-shock 
oscillations in the fluid. 

The standard SPH artificial viscosity contains two terms, one that is linear in the 
particle velocity differences, and another that is quadratic (Monaghan 1992). The linear 
term has user-specified coefficient a and tends to dominate for low Mach numbers, while 
the quadratic term has coefficient (3 and is important for higher Mach numbers. Typically, 
the values a ~ 1 and (3 ~ 2 are used and the shock front is spread over ~ 3 — 4 particle 
smoothing lengths (Hernquist & Katz 1989). This type of SPH viscosity can introduce 
shear into the flow, particularly through the linear term (Hernquist & Katz 1989; Monaghan 
1992). Since shear can affect the bar mode instability we want to keep the artificial viscosity 
coefficients, particularly a, as small as possible while still maintaining accuracy in the 
presence of the shock waves that occur as the ends of the bar and the spiral arms expand 
supersonically. We have chosen to use a = 0.25 and (3 = 1.0 as our standard values. 

TREESPH contains the option to use a version of the usual SPH artificial viscosity 
that reduces the amount of artificial viscosity in the presence of curl (Balsara 1989; Balsara, 
et al. 1989; Benz 1990). Tests carried out by Centrella & McMillan (1993) using TREESPH 
to calculate the gravitational radiation from the head-on collision of identical polytropes 
showed that the best results were obtained using this modified artificial viscosity. We have 
therefore chosen to use it here. 

In these simulations shocks occur in the outer regions of the model as the bar and 
spiral arms develop. The effects of changing the artificial viscosity coefficients can be seen 
by comparing models T3 and T4. Both of these runs have N = 15, 648 and started from 
the same initial equilibrium model; see Table |j. Run T3 has our standard values a = 0.25 
and (3 = 1.0, while T4 uses the larger values a = 1.0 and (3 = 2.0. Models T3 and T4 are 
very similar in visual appearance and in their bulk properties, having the same core values 
in the final state; see Table [5|. The behavior of the Fourier components of the density in 
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both models is also similar. Although the bar mode growth rate is slightly lower in T4, 
this may not be significant due to the element of subjectivity involved in calculating it. As 
expected, the energy conservation in T4 is improved due to the greater smoothing of shocks 
by the larger artificial viscosity. In addition, T4 requires significantly more CPU time than 
T3, due to the stability requirement on the particle timesteps in the presence of artificial 
viscosity (Hernquist & Katz 1989). 

The gravitational waveforms and luminosities for runs T3 and T4 are similar, except 
that the quantities in T4 have lower amplitudes than in T3. This can be compared with the 
results of Zhuge, Centrella, & McMillan (1994), who used TREESPH to simulate binary 
neutron star coalescence. In their models, the stars merge and coalesce into a rotating 
bar-like structure. Spiral arms form as mass is shed from the ends of the bar; the arms 
then expand and merge into a disk around the central object. They found that, during this 
spiral arm stage, the amplitude of the gravitational waveforms decreases as the amount of 
artificial viscosity is increased. 

For comparison, run T5 also has the same number of particles and began with the same 
initial state as T3, but has no explicit artificial viscosity, with a = (3 = 0. In this case, the 
kinetic energy of the particles passing through shocks is not correctly converted into heat, 
resulting in large post-shock oscillations. Table ^ shows that the energy conservation errors 
are more than twice as large as those for run T3. The growth of the Fourier components is 
similar, but the spiral arms spread out and merge more quickly after the bar mode peaks in 
T5 compared to T3. The final core w < R eq has less mass and angular momentum than in 
run T3. The values of the stability parameter r for both the core and the entire system are 
about the same for T5 and T3. However, T5 has a much larger kinetic energy due to the 
large post-shock oscillations. This gives T tota i/|W| = 0.31 for the whole system compared 
with T tot£ a/|W| = 0.26 for T3. The waveforms and luminosity profiles for model T5 are 
similiar to T3 until around the time that the amplitude of the bar mode reaches its peak 
value; afterwards, the profiles for T5 become more noisy. 

Finally, we compare models T6 and T3 to see the effects of starting with a different 
initial model. Run T6 uses the same values of a and (3 as T3 and N = 16, 000 particles; 
this is the model that was presented in Houser, Centrella, & Smith (1994). The initial 
conditions for this run were produced using the random or rejection method. As mentioned 
in § W\ this technique produces relatively large density fluctuations about the equilibrium 
solution. Thus the particles in the initial model are acted on by forces that can have large 
deviations from their expected values, which can lead to violent motions (Lucy 1977). 
Although TREESPH does perform some smoothing of the initial data (Hernquist & Katz 
1989), significant fluctuations remain. This results in energy conservation errors ~ 4.3%, 
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which are more than twice as large as those in T3. The amplitude of the bar mode in T6 
starts growing from a larger initial value and reaches a much broader peak before dropping 
to a local minimum again. In particular, there is a relatively short time interval during 
which we can define a linear growth region for In \C\ so we must use caution in interpreting 
the relatively large growth rate reported in Table [|. Run T6 also has more mass and 
angular momentum in the core at the final state than T3, although the stability parameter 
has the same value in both cases. Comparison of the gravitational wave data shows that the 
gravitational wave burst begins immediately in T6, compared to a starting time ~ IO^d for 
T3. Since T6 goes out to 26.9£d whereas T3 goes out to 35£d, the final states are roughly 
equivalent and can be compared meaningfully. The gravitational wave quantities for T6 
do show the signature of a rotating nonaxisymmetric core after the burst, although there 
is not as much detail as in T3. The amplitudes of the waveforms and the luminosities are 
both lower for T6. See Table H 

It is interesting to compare our results with those of DGTB, who evolved the dynamical 
bar instability using an SPH code with N = 2000 particles and smoothing lengths that are 
allowed to vary in time but are the same for all particles. Plots of the particle positions 
projected onto the x — y plane for the case r ~ 0.33 are visually similar to our Run Tl. 
They report that the system has r = 0.247 at the end of their run, which is essentially the 
same as our result for Tl. 

In summary, all the models run with TREESPH show the development of the bar and 
the spiral arm pattern. Models Tl - T4 and T7 conserve total energy to <^ 2%, with T5 and 
T6 conserving energy to ~ 4%. In all cases, angular momentum is conserved to ^ 0.1%. 
The bar mode growth rates for models T2 - T7 are the same to within ~ 6%; we attribute 
these differences largely to the element of subjectivity inherent in choosing the region of 
linear growth for <iln \ C\/dt. The anomalous growth rate found for Tl is due to the lack 
of a clearly defined linear growth region. The properties of the final cores w < R eq in the 
models are remarkably similar, especially if runs T5 (with no artificial viscosity) and T6 
(with a noisy initial state) are excluded. Finally, the amplitudes of the gravitational wave 
quantities increase as N increases; higher resolution runs, or models with non-equal-mass 
particles are needed to achieve convergence in these quantities. 



7. Comparison of Eulerian and SPH Results 

We turn now to a comparison between the results of the Eulerian and SPH codes. To 
accomplish this we focus on our standard models E2 and T7. There is no simple measure 
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for comparing the resolution of these two models since the underlying fluid descriptions are 
so different. In run E2, the fluid initially occupies ~ 14,400 zones; after flowing through 
the grid during the development of the bar mode, the fluid occupies ~ 66, 200 zones at the 
end of the run. In contrast, the fluid in run T7 is discretized into N = 32, 914 equal-mass 
particles. As the system evolves, these particles move through space, with the smoothing 
length of each particle continually adjusted to keep the number of nearest neighbors 
approximately constant. For our purposes we simply note that the number of particles used 
in T7 is comfortably between the initial and final number of zones occupied by the fluid in 
E2, and proceed with the comparison. 

Examination of Tables |IJ and |3] shows that T7 conserves both total energy and angular 
momentum better than does E2. Run T7 also uses a larger amount of CPU time. However, 
as demonstrated in § |5.2.|, the energy and angular momentum conservation for E2 both 



improve when relaxed vacuum conditions are used, although at the expense of a larger 
usage of CPU time. 

In both models, the growth of the m = 2 mode produces a rotating bar that develops 
spiral arms as mass is shed from the ends of the bar. The bar and spiral arms exert 
gravitational torques that cause angular momentum to be transported outward. The spiral 
arms expand supersonically and merge together, causing shock heating and dissipation in 
the disk surrounding the central core. The system remains highly flattened and evolves 
toward a nearly axisymmetric final state in both models. 

The density contours for E2 in Figure |l] and for T7 in Figure ^ can be compared 
directly, since the contours for T7 were made by interpolating the particle model onto the 
grid used for E2 and these figures use the same contour levels. The most striking visual 
difference lies in the fact that T7 expands much more than does E2. These differences in 
the amount of expansion account for the differences in the core [w < R eq ) masses and 
angular momenta given in Tables ||| and |5|. When the net transport of angular momentum 
is considered, the behavior of the models is more similar. At the final time, for example, in 
model E2 90% of the mass has 73% of the angular momentum and in T7 90% of the mass 
has 71% of the angular momentum. 

The growth of the bar mode is shown quantitatively in Figure [| for run E2 and in 
Figure ITS] for run T7. In both cases there is a relatively long period of linear growth for 



the bar mode amplitude ln|C|. The measured values in this linear growth region are 
dln\C\/dt = 0.58^ for E2 and d\n\C\/dt = 0.51^ for T7. The larger slope for E2 is 
closer to the analytic TVE value d\n\C\/dt = 0.728 ± 0.038^ given by TDM. As was 
discussed in § |5.2.| , numerical diffusion can cause the growth rate given by a simulation 
to be lower than the expected value (TDM). This suggests that the SPH code has more 
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numerical diffusion than the Eulerian code; further studies of this would be very useful. 
The eigenfrequency is a 2 = lSt^ 1 for E2 and a 2 = 1.9^ for T7; both of these values are 
within the range a 2 = 1.892 ± 0.094t5 1 given for the TVE result (TDM). 

In both runs, the bar mode amplitudes reach their maximum values, drop off to 
local minima, and then grow again. The bar mode in T7 reaches a higher value than in 
E2, but the peak is broader. In run E2 the amplitude reaches a second peak at a lower 
amplitude than the first, and then drops again. In both cases, the m = 4 mode starts 
growing after the bar mode, grows at a faster exponential rate, and then peaks at about the 
same time as the bar mode before dropping off; cf. Figures ||(c) and 0(c). For E2 we find 
d ]n(|C 4 |/|C |)M = and a 4 = 3.4^* and for T7 we get d\n(\C 4 \/\C \) /dt = l.lt D l 

and (T4 = In both runs we find that the pattern speeds for these two modes are 

about the same, indicating that the m = 4 mode is a harmonic of the m = 2 mode. 

One interesting difference between the models can be seen by comparing the behavior 
of the m — 1 and m = 3 Fourier components of the density. Figure |5] shows that these 



disturbances grow in run E2 whereas Figure 16 shows that they do not grow in run T7. We 



do not have an explanation for this behavior. However, given the importance of the m = 1 
mode in recent work (Bonnell 1994; Bonnell & Bate 1994; PDD), this question deserves 
further study 

The overall behavior of the stability parameter r and T tota i/|W| is similar in both E2 



and T7, as can be seen by comparing Figures |6| and [17]. In both cases, the final value of 
the stability parameter is in the range r s < r < We therefore expect that the system 
will develop a secular bar instability when the effects of gravitational radiation reaction are 
included in the hydrodynamical equations (cf. Lai & Shapiro 1995). 

The gravitational radiation in both models is dominated by a strong feature that 
corresponds to the initial growth of the bar mode, with the peak amplitudes of both the 
waveforms and luminosity being higher in E2 than in T7. The radiation emitted after 
the initial burst shows a secondary feature in both models. In E2 the radiation has a 
double-burst structure; this behavior is less distinct in T7, although the secondary feature 



is present. We showed in § |5.2.| that the gravitational radiation amplitudes for E2 decrease 
somewhat if we halve the number of angular zones or change the treatment of the vacuum 
zones to allow the model to expand more. Also, § |6.2.| showed that the gravitational wave 
amplitudes for the SPH models increase and approach the E2 amplitude as the particle 
number increases, and that they have not yet converged for this set of runs. Further work 
with both codes, including higher resolution runs, is needed to resolve these issues 
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8. Conclusions 

We have carried out 3-D numerical simulations of the dynamical bar instability in 
a rapidly rotating star and the resulting gravitational radiation using both an Eulerian 
finite-difference code with a cylindrical grid and monotonic advection, and an SPH code with 
variable smoothing lengths and a hierarchical tree method for calculating the gravitational 
acceleration. The star is initially modeled as a polytrope with index n = 3/2 and r m 0.30. 
In both codes the gravitational field is purely Newtonian and the gravitational radiation is 
calculated using the quadrupole formula. The back reaction of the gravitational radiation 
on the fluid is not included. 

In both codes the dynamical instability of the m = 2 mode produces a rotating bar-like 
structure. Spiral arms develop as mass is shed from the ends of the bar, and gravitational 
torques cause angular momentum to be transported outward. The spiral arms expand 
supersonically and merge, causing shock heating in the outer regions. At the end of the 
simulations, both codes agree that the system consists of a nearly axisymmetric central core 
surrounded by an extended disk, and has r ~ 0.25. 

It is interesting to compare our results with those from the earlier study by DGTB 
that considered rapidly rotating n = 3/2 polytropes with r ps 0.33 and r 0.38. The 
Eulerian codes used in that work were quite diffusive. They also did not solve an energy 
equation and thus had no means of handling self-consistently the shocks that form in the 
outer regions. Also, a certain amount of mass was allowed to leave the grid. In our Eulerian 
runs, shocks are handled using an artificial viscosity and no mass is allowed to leave the 
grid. Our main difficulty in the Eulerian code is with the expansion of the model into the 
vacuum. We believe this can be alleviated with the use of better vacuum conditions. 

The SPH code used by DGTB had a smoothing length that was the same for all 
particles. Their runs were also limited to a fairly small number of particles. The use of 
variable smoothing lengths and the hierarchical tree method for calculating the gravitational 
accelerations in TREESPH allows us to evolve more particles. All the SPH models have no 
grid and expand freely. However, the fluid description can break down in the low density 
outer regions due to lack of resolution. 

The simulations of DGTB showed the development of the bar instability with spiral 
arms and transport of angular momentum. However, their codes differed in the final 
outcome of the models in the low density outer regions. In their longest Eulerian run, most 
of the low density material formed a fairly narrow ring around the central remnant. (This 
outcome was also seen by Williams & Tohline (1988) for polytropes having n = 0.8 and 
n = 1.8 with r = 0.31.) However, in their SPH runs the low density material formed an 
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extended disk. 

We are not sure what causes these differences in the final outcome. Further study, 
perhaps including longer runs with our codes, is needed to resolve these issues. 

Interestingly, DGTB concluded that their SPH runs had less numerical diffusion than 
their Eulerian simulations. We believe this is due to the low-order differencing in their 
Eulerian codes. In our study, the use of higher-order differencing and monotonic advection 
in the Eulerian code resulted in much less numerical diffusion than seen in the earlier 
studies. (See also PDD.) Comparison of the bar growth rates for our standard Eulerian and 
SPH runs suggests that the SPH code has more numerical diffusion. 

We agree with DGTB that the SPH code is easier to implement and use. We find that 
the bulk properties of the model can be obtained at very low cost using iV ~ 2000 particles, 
although more particles are needed for a good measure of the bar mode growth rate. We 
did not carry out any very low resolution studies with our Eulerian code. However, we find 
that the peak amplitudes of the gravitational radiation quantities in both the Eulerian and 
SPH codes increase as the resolution is increased. Overall, for comparable resolution, the 
cost of the Eulerian and SPH runs is similar. 

We note that although the version of SPH used here allows the particle smoothing 
lengths to vary in both space and time, the terms describing these changing smoothing 
lengths are not explicitly incorporated into the dynamical equations (Hernquist & Katz 
1989). Although this situation is typical of most SPH codes currently used in astrophysics, 
it does present a potentially serious deficiency in the method. Recent work to incorporate 
these terms into the dynamical equations (e.g. Nelson & Papaloizou 1995) may lead to 
substantial improvements in the SPH method. 

Finally, the suitability of a particular numerical method must be determined in the 
context of the astrophysical system being modeled, as well as the resources available to the 
investigators. We hope that this study, along with related work comparing the results of 
Eulerian and SPH hydrodynamic codes in modeling stellar collisions (Davies, et al. 1993) 
and the growth of structure in cosmology (Kang, et al. 1994), will help others to find the 
best approach for their applications. 

We thank S. Clancy, S. Cranmer and S. McMillan for helpful conversations, and L. 
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comments of the referee R. Durisen, which contributed to improving this paper. This work 
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the Cray C90 at Pittsburgh Supercomputing Center under grant PHY910018P. 
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Run 




N z 


N v 


-^star,i 


N s tax,f 


Ti 


time 


Ei 




Ji-jf 
J\ 


CPU 


El 


64 


32 


32 


7200 


24,600 


0.30 


34t D 


.035 


.072 


8.2 hr 


E2 


64 


32 


64 


14,400 


66,200 


0.30 


34t D 


.043 


.086 


29.7 hr 


E2' 










41,500 




23.9t D 


.020 


.049 


20.5 hr 


E3 


64 


32 


64 


14,400 


110,600 


0.30 


23.9£ D 


.015 


.034 


29.5 hr 



Table 1: Properties of the Eulerian models. iV ro , N ZJ and N v are the number of grid zones in 
the zu, z, and ip directions, respectively. iV s tar is the approximate number of zones occupied 
by the matter distribution. The subscripts "i" and "f" denote the initial and final states 
of the model. The stability parameter at the initial time, calculated on the Eulerian grid, 
is T{. The duration of the run is measured in units of the dynamical time in and listed in 
the column labeled "time". All models were run on a Cray C90; the amount of CPU time 
used is given for the duration of the run. The quantities given in the row labeled E2' are the 
values for model E2 at the intermediate time 23.9£d- Model E3 is the same as E2 except for 
the use of relaxed constraints in the vacuum conditions. 



Run 


d\n\C\/dt 


0"2 


-^core.f 


"^core,f 


7"core,f 


Tf 


El 


0.53 


1.7 


94% 


78% 


0.18 


0.19 


E2 


0.58 


1.8 


96% 


86% 


0.24 


0.24 


E2' 






96% 


87% 


0.25 


0.26 


E3 


0.57 


1.8 


95% 


83% 


0.25 


0.26 



Table 2: Hydrodynamical and bar mode results for the Eulerian models. The quantities in 
the row labeled E2' are the values for run E2 at the intermediate time 23.9to; cf. Table [l]. The 
bar growth rate din \C\/dt and the eigenfrequency o"2 are calculated for the ring w = 0.362i? eq 
in zone j = 10 in the equatorial plane; both are in units of . The core refers to matter 
within cylindrical radius w = R eq , where R eq is the equatorial radius of the initial model, 
and the subscript "f" denotes the final state of the model. Note that J CO re,f is normalized 
using the value of the total angular momentum of the system at the final time, which is less 
than the initial value due to non-conservation; see Table |l|. 
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Run max \rh\ max L [AE/M){ ma,xdJ z /dt (AJ/Jo)f 

El 0.63 0.20 0.68 0.18 1.6 

E2 0.68 0.21 0.93 0.20 2.3 

E2' 0.77 2.0 

E3 0.62 0.17 0.63 0.16 1.6 



Table 3: Gravitational wave results for the Eulerian models. We use c = G = 1. The 
quantities in the row labeled E2' are the values for run E2 at the intermediate time 23.9£d; 
cf. Table |l|. The values of the peak gravitational wave amplitudes \rh\ are in units of M 2 /R eq . 
The values of the maximum luminosity L are in units of (M/R eq ) 5 , and the values of the 
total energy emitted during the duration of the run (A£ , /M) f are in units of [Mj R eq ) 7 ^ 2 . 
The values of the maximum dJ z /dt are in units of M(M/R cq ) 7 ^ 2 . The quantity (AJ Z /J ){ 
is the total angular momentum emitted as gravitational radiation normalized by the initial 
total angular momentum, and has units of (M/i? eq ) 5//2 . 



Run 


N 


type 


a 


(3 




n 


time 


h 


-E { 
Ei 




./ 


-J{ 
Ji 


CPU 


Tl 


2061 


cold 


0.25 


1.0 





305 


35t D 





022 




.001 


1.01 hr 


T2 


8728 


cold 


0.25 


1.0 





308 


35t D 





018 




.001 


8.65 hr 


T3 


15,648 


cold 


0.25 


1.0 





314 


35i D 





018 




.001 


16.9 hr 


T4 


15,648 


cold 


1.0 


2.0 





314 


35t D 





011 




.001 


23.0 hr 


T5 


15,648 


cold 











314 


35i D 





043 




.001 


16.7 hr 


T6 


16,000 


random 


0.25 


1.0 





299 


26.9t D 





043 




.001 


11.6 hr 


T7 


32,914 


cold 


0.25 


1.0 





316 


35t D 





018 




.001 


44.6 hr 



Table 4: Properties of the SPH models. All models used the "cold" initial conditions except 
T6; the initial conditions for this model were produced using the random method. N is the 
number of particles in the model, and a and f3 are, respectively, the coefficients of the linear 
and quadratic terms in the artificial viscosity. The other quantities are defined as in Table |I[ 
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Run 


d\n\C\/dt 


(T 2 


-^core,f 


J core,f 


Tcore,f 


Tf 


1 1 


n p.p. 
u.oo 


Z. 1 


yi /O 


to /o 


n OA 


u.zo 


T2 


0.51 


1.9 


92% 


74% 


0.24 


0.26 


T3 


0.53 


1.9 


92% 


73% 


0.25 


0.26 


T4 


0.51 


1.9 


92% 


73% 


0.25 


0.26 


T5 


0.54 


1.9 


89% 


67% 


0.24 


0.26 


T6 


0.54 


1.8 


95% 


81% 


0.25 


0.26 


T7 


0.51 


1.9 


90% 


72% 


0.25 


0.26 



Table 5: Hydrodynamical and bar mode results for the SPH models. The quantities are 
defined as in Table The large growth rate for Tl is anomalous; see the text for details. 



Run 


max \rh\ 


max L 


(AE/M)f 


max dJ z /dt 


(AJ z /J ) f 


Tl 


0.26 


0.030 


0.21 


0.015 


0.24 


T2 


0.43 


0.078 


0.46 


0.039 


0.54 


T3 


0.47 


0.091 


0.55 


0.045 


0.66 


T4 


0.41 


0.066 


0.39 


0.035 


0.48 


T5 


0.48 


0.10 


0.52 


0.050 


0.59 


T6 


0.38 


0.051 


0.31 


0.027 


0.42 


T7 


0.53 


0.12 


0.87 


0.060 


1.0 



Table 6: Gravitational wave results for the SPH models. The quantities are defined as in 
Table |3|. 



-33- 



REFERENCES 

Abramovici, A., et al. 1992, Science, 256, 325 
Balsara, D. 1989, Ph.D. Thesis, University of Illinois 

Balsara, D., Norman, M., Bicknell, G., & Gingold, R. 1989, unpublished manuscript 
Barnes, J. & Hut, P. 1986, Nature, 324, 446 

Benz, W. 1990, in Numerical Modeling in Stellar Pulsation: Problems and Prospects, ed. 
J. Buchler (Dordrecht: Kluwer) 

Bodenheimer, P. & Ostriker, J. 1973, ApJ, 180, 159 

Bonnell, I. 1994, MNRAS, 269, 837 

Bonnell, I. & Bate, M. 1994, MNRAS, 271, 999 

Bowers, R. & Wilson, J. R. 1991, Numerical Modeling in Applied Physics and Astrophysics 
(Boston: Jones & Bartlett) 

Bradaschia, C, et al. 1990, Nucl. Instrum. Methods Phys. Research, A289, 518 

Centrella, J. & McMillan, S. 1993, ApJ, 416, 719 

Chandrasekhar, S. 1969, Ellipsoidal Figures of Equilibrium (New Haven: Yale Univ. Press) 

Clancy, S. P. 1989, Ph.D. Thesis, University of Texas, Austin 

Colella, P. & Woodward, P. 1984, J. Comp. Phys., 82, 64 

Davies, M., Ruffert, M., Benz, W., & Miiller, E. 1993, A&A, 272, 430 

Durisen, R., Gingold, R., Tohline, J., & Boss, A. 1986, ApJ, 305, 281 (DGTB) 

Durisen, R. & Tohline, J. 1985, in Protostars & Planets II, eds. D. Black & M. Matthews 
(Tucson: University of Arizona Press) 

Finn, L. S. & Evans, C. R. 1990, ApJ, 351, 588 

Gingold, R. & Monaghan, J. 1977, MNRAS, 181, 375 

Hachisu, I. 1986, ApJS, 61, 479 

Hawley, J., Smarr, L., & Wilson, J. 1984, ApJS, 55, 211 
Hernquist, L. 1987, ApJS, 64, 715 
Hernquist, L. 1993, ApJ, 404, 717 
Hernquist, L. & Katz, N. 1989, ApJS, 70, 419 
Houser, J. & Centrella, J. 1995, in preparation 



-34- 



Houser, J., Centrella, J., & Smith, S. 1994, Phys.Rev.Lett, 72, 1314 

Imamura, J., Friedman, J., & Durisen, R. 1985, ApJ, 294, 474 

Imamura, J., Toman, J., Durisen, R., Pickett, B., & Yang, S. 1995, ApJ, 444, 363 

Kang, H., et al. 1994, ApJ, 430, 83 

Lai, D. & Shapiro, S. 1995, ApJ, 442, 259 

Landau, L. D. & Lifshitz, E. M. 1959, Fluid Mechanics (Reading, MA: Addison Wesley) 
Lucy, L. 1977, AJ, 82, 1013 
Managan, R. 1985, ApJ, 294, 463 

Meijerink, J. A. k Van Der Vorst, H. A. 1981, J. Comp. Phys., 44, 134 

Misner, C, Thorne, K., & Wheeler, J. 1973, Gravitation (New York: Freeman) 

Monaghan, J. 1992, ARA&A, 30, 543 

Monaghan, J. & Lattanzio, J. 1985, A&A, 149, 135 

Nelson, R. & Papaloizou, J. 1995, MNRAS, in press 

Norman, M. L., Wilson, J. R., & Barton, R. 1980, ApJ, 239, 968 

Norman, M. L. & Winkler, K.-H., 1986, in Astrophysical Radiation Hydrodynamics, eds. 
K.-H. Winkler & M. L. Norman (Dordrecht: Reidel) 

Ostriker, J. & Bodenheimer, P. 1973, ApJ, 180, 171 

Ostriker, J. & Mark, J. 1968, ApJ, 151, 1075 

Pickett, B., Durisen, R., & Davis, G. 1996, ApJ, in press (PDD) 

Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. 1992, Numerical 
Recipes, 2 nd Edition (New York: Cambridge University Press) 

Schutz, B. 1986, in Dynamical Spacetimes and Numerical Relativity, ed. J. Centrella (New 
York: Cambridge University Press) 

Schutz, B. 1989, Class. Quantum Gravity 6, 1761 

Shapiro, S. & Teukolsky, S. 1983, Black Holes, White Dwarfs, and Neutron Stars (New 
York: Wiley) 

Smith, S. C. 1993, Ph.D. Thesis, Drexel University 

Smith, S. C. & Centrella, J. M. 1992, in Approaches to Numerical Relativity, ed. R. 
d'Inverno (New York: Cambridge University Press) 

Smith, S. C, Centrella, J. M., & Clancy, S. P. 1994, ApJS, 94, 789 



-35 - 



Tassoul, J.-L. 1978, Theory of Rotating Stars (Princeton: Princeton University Press) 
Tassoul, J.-L. & Ostriker, J. P. 1968, ApJ, 154, 613 

Thorne, K. 1995, in Relativistic Cosmology, Proceedings of the Eighth Nishinomiya- Yukawa 
Memorial Symposium, ed. M. Sasaki (Tokyo: Universal Academy Press) 

Tohline, J., Durisen, R., & McCollough, M. 1985, ApJ, 298, 220 (TDM) 

Williams, H. & Tohline, J. 1987, ApJ, 315, 594 

Williams, H. & Tohline, J. 1988, ApJ, 334, 449 

Wilson, J. R. 1979, in Sources of Gravitational Radiation, ed. L. Smarr (Cambridge: 
Cambridge University Press) 

Zhuge, X., Centrella, J., & McMillan, S. 1994, Phys.Rev.D, 50, 6247. 



This preprint was prepared with the AAS IATeX macros v3.0. 



-36- 



Fig. 1. — Density contours in the equatorial plane for model E2. The contour levels are 
the same in all frames. The contours are spaced a factor of 10 apart, and go down to 3 
decades below the maximum (central) density at the initial time t = 0. Since the maximum 
density increases slightly during the run, there are 4 contours shown for the later times, 
with the innermost contour being at the initial maximum density. The model rotates in the 
counterclockwise direction. The circular boundary of the plotted region is set at w — 2.5R eq . 

Fig. 2. — The mass fraction m(w) is shown for model E2 at the initial time (dot-dashed 
line), the intermediate time t = 21.1£d (dashed line), and the final time t = 34£d (solid line). 
The total mass is M. 

Fig. 3. — The angular momentum J{w) is shown for the model E2 at the initial time (dot- 
dashed line), the intermediate time t = 21. lt^ (dashed line), and the final time t = 34tc 
(solid line). Note that J(w) is normalized by the total angular momentum in the system at 
the time. 

Fig. 4. — The growth of the bar mode for model E2. The bar amplitude In |C| is shown 
versus time for the ring w = 0.362i? cq in zone j = 10 in the equatorial plane z = 0. The 
growth rate in the region where In \C\ is growing linearly with time is <iln \C\/dt = 0.58^. 

Fig. 5. — The growth of various Fourier components of the density for model E2. The 
amplitudes have been calculated for the same ring used in Figure f|, and are shown versus 
time. In each case, the bar mode amplitude is plotted as a solid line for comparison, (a) 
m = 1 (b) m = 3 (c) m = 4 

Fig. 6. — The behavior of T/|W| is shown as a function of time for model E2. The solid line 
shows the stability parameter r = T rot /|W| and the dashed line shows T to tai/|W|- 

Fig. 7. — Gravitational waveforms for an observer on the axis at 9 = and ip = for model 
E2. We use c = G = 1. (a) rh + (b) rh x 
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Fig. 8. — Various gravitational wave quantities are shown for run E2. We use c = G = 1. (a) 
Gravitational wave luminosity L. (b) The energy AE/M emitted as gravitational radiation, 
(c) The rate dJ z /dt at which angular momentum is carried away by the waves, (d) The 
angular momentum A J/ J lost to gravitational radiation. Here, Jo is the initial total angular 
momentum. 

Fig. 9. — Density contours in the equatorial plane at t — 23.9£d- The contour levels are the 
same as in Figure 0. (a) Model E2 (b) Model E3 (with relaxed vacuum conditions) 

Fig. 10. — Particle positions are shown projected onto the equatorial plane for various times 
in the evolution of model T7. The vertical axis is y/R eq and the horizontal axis is x/R eq . 
The system rotates in the counterclockwise direction. 

Fig. 11. — Particle positions are shown projected onto the x — z plane at the final time 
t = 35£d of model T7. Shock heating during the spiral arm expansion has caused the disk 
to puff up. 

Fig. 12. — Density contours in the equatorial plane for model T7. The density has been 
interpolated onto the cylindrical grid used for model E2, and the contour levels are the same 
as those used in Figure [IJ. 

Fig. 13. — The mass fraction m{w) is shown for model T7 at the initial time (dot-dashed 
line), the intermediate time t = 18.2£d (dashed line), and the final time t = 35£d (solid line). 
Only the region w < 5R cq is plotted. The total mass is M. 

Fig. 14. — The angular momentum J(w) is shown for the model T7 at the initial time (dot- 
dashed line), intermediate time t = 18.2t D (dashed line), and the final time t = 35£d (solid 
line). Only the region w < 5R eq is plotted. The total initial angular momentum is Jo- 

Fig. 15. — The growth of the bar mode for model T7. The particle model is interpolated onto 
the cylindrical grid used for model E2, and the bar amplitude is calculated as in Figure |j. The 
growth rate in the region where In \C\ is growing linearly with time is din \C\/dt = O.Slt^ 1 - 

Fig. 16. — The growth of various Fourier components of the density for model T7. In each 
case, the bar mode amplitude is plotted as a solid line for comparison, (a) m = 1 (b) m = 3 
(c) m — 4 
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Fig. 17. — The behavior of is shown as a function of time for model T7. The solid 

line shows the stability parameter r = T rot /|W| and the dashed line shows T total /|VF|. The 
initial oscillations in T/|W| are caused by adjustments in the model due to residual noise in 
the initial conditions. 



Fig. 18. — Gravitational waveforms for an observer on the axis at 6 = and (f = for model 
T7. We use c = G = 1. (a) rh+ (b) rh x 



Fig. 19. — Various gravitational wave quantities are shown for run T7. We use c — G — 1. (a) 
Gravitational wave luminosity L. (b) The energy AE/M emitted as gravitational radiation, 
(c) The rate dJ z /dt at which angular momentum is carried away by the waves, (d) The 
angular momentum A J / J lost to gravitational radiation. 



